2.1.1 공기질 데이터 불러오기

구득한 엑셀데이터를 ’readxl’라이브러리를 사용하여 불러온 후, 측정소명이 ’로’로 끝나는 행만 가져왔습니다.

AQraw <- read_xlsx("./data/일별평균대기오염도_2018.xlsx")
AQroad <- subset(AQraw, grepl(glob2rx("*로"), 측정소명))

head(AQroad)

# A tibble: 6 x 8
  측정일시 측정소명 `이산화질소농도(ppm)`~ `오존농도(ppm)` `이산화탄소농도(ppm)`~
     <dbl> <chr>               <dbl>           <dbl>            <dbl>
1 20180101 강남대로            0.04            0.007              0.8
2 20180101 강변북로            0.033           0.008              0.6
3 20180101 공항대로            0.045           0.007              0.6
4 20180101 도산대로            0.037           0.009              0.6
5 20180101 동작대로            0.051           0.005              0.7
6 20180101 신촌로              0.045           0.005              1  
# ... with 3 more variables: `아황산가스(ppm)` <dbl>, `미세먼지(㎍/㎥)` <dbl>,
#   `초미세먼지(㎍/㎥)` <dbl>

2.1.2 관측지점 자료와의 조인

구득한 관측지점을 시각화하고자 관측지점 주소 데이터를 조인하였습니다.

AQlocation <- read_xls("./data/station_list.xls", sheet = 1)[-c(1, 2), -c(1, 4, 5)]
AQlocation[5, 1] <- "동작대로"
names(AQlocation) <- c("측정소명", "주소")

head(AQlocation)

# A tibble: 6 x 2
  측정소명 주소                                                            
  <chr>    <chr>                                                           
1 강남대로 서울특별시 서초구 강남대로 201 서초구민회관 앞 중앙차로 (양재동)
2 강변북로 서울 성동구 강변북로 257 한강사업본부 옆                        
3 공항대로 서울 강서구 마곡동 727-1091 마곡역 중앙차로정류장 옆            
4 도산대로 서울 강남구 도산대로 104 신사역2번출구 앞                       
5 동작대로 서울 동작구 동작대로 144 이수역 북단 버스중앙차로               
6 시흥대로 서울 금천구 독산동 996-9 한양수자인아파트 앞                    

locate_AQroad <- inner_join(AQroad, AQlocation, by = "측정소명")
locate_AQroad

# A tibble: 4,362 x 9
   측정일시 측정소명 `이산화질소농도(ppm)`~ `오존농도(ppm)` `이산화탄소농도(ppm)`~
      <dbl> <chr>               <dbl>           <dbl>            <dbl>
 1 20180101 강남대로            0.04            0.007              0.8
 2 20180101 강변북로            0.033           0.008              0.6
 3 20180101 공항대로            0.045           0.007              0.6
 4 20180101 도산대로            0.037           0.009              0.6
 5 20180101 동작대로            0.051           0.005              0.7
 6 20180101 신촌로              0.045           0.005              1  
 7 20180101 영등포로            0.049           0.007              0.6
 8 20180101 정릉로              0.05            0.007              0.8
 9 20180101 종로               NA              NA                 NA  
10 20180101 천호대로            0.045           0.007              0.6
# ... with 4,352 more rows, and 4 more variables: `아황산가스(ppm)` <dbl>,
#   `미세먼지(㎍/㎥)` <dbl>, `초미세먼지(㎍/㎥)` <dbl>, 주소 <chr>

2.1.3 조인을 위한 인덱스 생성 1

같은 위치상에 존재하는 대기질 관측소와 교통량 정보 수집 장치의 정보를 구할 수 없어, 비슷한 위치상의 것들로 대체하였습니다. 주소가 달라 조인 인덱스로 사용할 수 없었기에, 날짜와 도로 명을 합친 인덱싱을 위한 하나의 새로운 열을 생성하였습니다.

index <- vector(mode = "character", length(locate_AQroad$측정일시))
for (h in 1:length(locate_AQroad$측정일시)) {
    index[h] <- paste(locate_AQroad$측정일시[h], locate_AQroad$측정소명[h])
}
AQ <- cbind(locate_AQroad, index)

datatable(AQ, caption = "공기질 데이터", filter = "top")

2.1.4 교통량 데이터 불러오기

구득한 데이터는 월별로 되어있어, for문을 사용하여 불러온 뒤 열을 기준으로 마지막 행 뒤에 차례로 합하였습니다.

fi <- paste0("./data/", list.files(path = "./data", pattern = "* 조사자료.xlsx"))
head(fi, n = 3L)

[1] "./data/01월 서울시 교통량 조사자료.xlsx"
[2] "./data/02월 서울시 교통량 조사자료.xlsx"
[3] "./data/03월 서울시 교통량 조사자료.xlsx"

for (g in 1:12) {
    assign(paste0("T", g, "raw"), read_xlsx(fi[g]))
}


Traw <- bind_rows(T1raw, T2raw, T3raw, T4raw, T5raw, T6raw, T7raw, T8raw, T9raw, 
    T10raw, T11raw, T12raw)

Traw

# A tibble: 105,304 x 30
     일자 요일  지점명 지점번호 구분  방향  `0시` `1시` `2시` `3시` `4시` `5시`
    <dbl> <chr> <chr>  <chr>    <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2.02e7 월    성산로(금~ A-01     유입  봉원고가~   657   617   405   279   266   365
 2 2.02e7 화    성산로(금~ A-01     유입  봉원고가~   382   256   166   171   266   751
 3 2.02e7 수    성산로(금~ A-01     유입  봉원고가~   641   492   310   238   329   691
 4 2.02e7 목    성산로(금~ A-01     유입  봉원고가~   661   483   357   279   357   644
 5 2.02e7 금    성산로(금~ A-01     유입  봉원고가~   764   524   361   282   339   672
 6 2.02e7 토    성산로(금~ A-01     유입  봉원고가~   807   594   486   412   372   521
 7 2.02e7 일    성산로(금~ A-01     유입  봉원고가~   617   478   354   264   230   316
 8 2.02e7 월    성산로(금~ A-01     유입  봉원고가~   448   278   231   205   256   715
 9 2.02e7 화    성산로(금~ A-01     유입  봉원고가~   641   467   304   256   312   679
10 2.02e7 수    성산로(금~ A-01     유입  봉원고가~   589   451   291   262   295   569
# ... with 105,294 more rows, and 18 more variables: `6시` <dbl>, `7시` <dbl>,
#   `8시` <dbl>, `9시` <dbl>, `10시` <dbl>, `11시` <dbl>, `12시` <dbl>,
#   `13시` <dbl>, `14시` <dbl>, `15시` <dbl>, `16시` <dbl>, `17시` <dbl>,
#   `18시` <dbl>, `19시` <dbl>, `20시` <dbl>, `21시` <dbl>, `22시` <dbl>,
#   `23시` <dbl>

2.1.5 두 데이터간 상응지점 지정

위에서도 언급하였듯이, 같은 위치상의 대기질 관측소와 교통량 수집기 데이터를 구득하지 못하였기에, 명칭과 주소를 대조하여 아래와 같이 11개의 서로 비슷한 지점을 찾을 수 있었습니다.

해당 작업은 if문을 사용하여 교통량 데이터의 지점에 상응하는 지점번호(코드)를 기준으로 대조를 통해 찾은 관측소들을 교통량 데이터 데이터프레임에 추가하였습니다.

# make specific correspondence w. AQ dataframe

a <- vector(mode = "character", length(Traw$지점번호))

for (i in 1:length(Traw$지점번호)) {
    if (Traw$지점번호[i] == "A-13") {
        # 1
        a[i] <- "한강대로"
    } else if (Traw$지점번호[i] == "A-15") {
        # 2
        a[i] <- "종로"
    } else if (Traw$지점번호[i] == "D-09") {
        # 3
        a[i] <- "정릉로"
    } else if (Traw$지점번호[i] == "D-10") {
        # 4
        a[i] <- "화랑로"
    } else if (Traw$지점번호[i] == "D-21") {
        # 5
        a[i] <- "공항대로"
    } else if (Traw$지점번호[i] == "D-24") {
        # 6
        a[i] <- "시흥대로"
    } else if (Traw$지점번호[i] == "D-25") {
        # 7
        a[i] <- "영등포로"
    } else if (Traw$지점번호[i] == "D-31") {
        # 8
        a[i] <- "동작대로"
    } else if (Traw$지점번호[i] == "D-35") {
        # 9
        a[i] <- "강남대로"
    } else if (Traw$지점번호[i] == "D-43") {
        # 10
        a[i] <- "도산대로"
    } else if (Traw$지점번호[i] == "F-02") {
        # 11
        a[i] <- "강변북로"
    } else {
        # else
        a[i] <- "NULL"
    }
}

cat_traffic <- cbind(Traw, a)

# remove null value in column a
cat1_traffic <- cat_traffic[!(cat_traffic$a == "NULL"), ]

2.1.6 일별 전체 합 구하기

head(cat1_traffic[, c(7:30)])

    0시 1시 2시 3시 4시 5시  6시  7시  8시  9시 10시 11시 12시 13시 14시 15시
745 514 458 258 288 194 360  549  535  783  883 1053 1218 1360 1383 1350 1411
746 326 287 189 176 230 539 1230 1976 2223 2164 1877 1936 1714 2067 1899 1911
747 565 406 300 211 223 519 1100 1986 2190 2124 2076 2101 2026 2037 1989 2068
    16시 17시 18시 19시 20시 21시 22시 23시
745 1287 1192  986  881  809  777  727  575
746 1835 1815 1809 1504 1219 1092  996  761
747 1842 1742 1729 1638 1193 1091  953  756
 [ reached 'max' / getOption("max.print") -- omitted 3 rows ]

# order by date
order_traffic <- cat1_traffic[c(order(cat1_traffic$일자)), ]

교통량 데이터가 하루 단위가 아닌 시간 단위로 되어 있고, 유출입이 나뉘어 있어, 하루 단위인 대기질 데이터와 맞춰줄 필요가 있었습니다.

위에서 먼저 시간순으로 정렬하여 지점별 유출입 행끼리 마주보게 하였습니다.

다음은 반복문을 사용하여 일별 총합을 구하는 과정입니다.

dsum <- vector(mode = "numeric", length(order_traffic$a))
for (j in 1:length(order_traffic$a)) {
    tmp <- vector(mode = "numeric", 24)
    for (k in 7:30) {
        tmp[k - 6] <- order_traffic[j, k]
    }
    dsum[j] <- sum(tmp)
}
iosum_traffic <- cbind(order_traffic, dsum)

head(iosum_traffic)

        일자 요일           지점명 지점번호 구분           방향 0시 1시 2시 3시
745 20180101   월 세종대로(서울역)     A-13 유입 서울역->숭례문 514 458 258 288
776 20180101   월 세종대로(서울역)     A-13 유출 숭례문->서울역 806 632 345 256
    4시 5시 6시 7시 8시 9시 10시 11시 12시 13시 14시 15시 16시 17시 18시 19시
745 194 360 549 535 783 883 1053 1218 1360 1383 1350 1411 1287 1192  986  881
776 257 367 461 454 691 980  972  968 1179 1251 1230 1328 1388 1334 1172 1053
    20시 21시 22시 23시        a  dsum
745  809  777  727  575 한강대로 19831
776 1022  944  716  608 한강대로 20414
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]

위처럼 지점별로 유출입 행이 마주보는 데이터 프레임을 가지고, 일일합 열을 2개씩 합하는 반복문을 사용하여 합하였습니다.

# sum dsum(io)
sumdsum <- vector(mode = "numeric", length(iosum_traffic$dsum)/2)
for (l in 1:as.numeric(length(iosum_traffic$dsum)/2)) {
    if (l == 1) {
        sumdsum[l] <- iosum_traffic$dsum[l] + iosum_traffic$dsum[l + 1]
    } else {
        sumdsum[l] <- iosum_traffic$dsum[l * 2 - 1] + iosum_traffic$dsum[l * 2]
    }
}

gyotong <- iosum_traffic[c(seq(from = 1, to = length(iosum_traffic$구분), by = 2)), 
    ]
traffic_final <- cbind(gyotong, sumdsum)[, -c(2, 4, 5, 6, 32)]

head(traffic_final)

        일자           지점명 0시 1시 2시 3시 4시 5시 6시 7시 8시 9시 10시 11시
745 20180101 세종대로(서울역) 514 458 258 288 194 360 549 535 783 883 1053 1218
869 20180101  종로(종로3가역)  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA   NA   NA
    12시 13시 14시 15시 16시 17시 18시 19시 20시 21시 22시 23시        a
745 1360 1383 1350 1411 1287 1192  986  881  809  777  727  575 한강대로
869   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA     종로
    sumdsum
745   40245
869      NA
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]

2.1.7 조인을 위한 인덱스 생성 2

대기질 데이터와 마찬가지로 날짜와 지정한 상응지점 값을 합하여 인덱싱을 위한 새로운 열을 생성하였습니다.

# generate index column
index <- vector(mode = "character", length(traffic_final$지점명))
for (m in 1:length(traffic_final$지점명)) {
    index[m] <- paste(traffic_final$일자[m], traffic_final$a[m])
}
index_traffic <- cbind(traffic_final, index)

datatable(index_traffic, caption = "서울시 일별 교통량", filter = "top")

2.1.8 대기질-교통량 데이터 조인

생성한 인덱스 열을 바탕으로 두 데이터를 조인하였습니다. inner-join을 사용하여 두 인덱스 열에 존재하지 않는 데이터는 자연스럽게 사라지게 하였습니다. 조인 후에 필요하지 않을 것이라 판단되는 열은 제거하였습니다.

# join
join <- inner_join(AQ, index_traffic, by = "index")[, c(1:9, 12, 38)]

2.1.9 결측치 제거

데이터의 결측치를 제거하였습니다

# remove missing value row
data_refined <- filter(join, !is.na(sumdsum) & !is.na(`미세먼지(㎍/㎥)`))

datatable(data_refined, caption = "서울시 일별 교통량", filter = "top")

2.2.1 지형자료로 변환을 위한 손질

지형자료 속성 테이블로의 변환을 위한 전처리 과정입니다.

# defining specific time to extract each observation station from time series
# data.
AQinfo <- subset(data_refined, 측정일시 == 20181231)

# traffic collector locations
trafinfo <- read_xlsx("./data/서울시 교통량 수집정보.xlsx", sheet = 2)[-c(136, 137), 
    -c(3, 4)] %>% cbind(read_excel("./data/서울시 교통량 수집정보.xlsx", sheet = 2, 
    col_types = "numeric")[-c(136, 137), c(3, 4)]) %>% rename(지점명 = "지점명칭")
trafinfo1 <- inner_join(trafinfo, AQinfo, by = "지점명")[, c(1, 2, 8, 7)]

2.2.2 지오코딩

R에서 지오코딩을 하기 위하여 Google의 API를 사용하였습니다. (API 키는 개인 키라 가렸습니다.)

register_google("구글API키주세요!")  # 과제 게시시 해당 청크 삭제할 것!
AQinfo1 <- as.data.frame(AQinfo)  # tibble will not work

gc <- mutate_geocode(AQinfo1, location = 주소, source = "google")

gc

  측정일시 측정소명 이산화질소농도(ppm) 오존농도(ppm) 이산화탄소농도(ppm)
1 20181231 강남대로               0.063         0.004                 1.2
2 20181231 강변북로               0.057         0.004                 0.6
3 20181231 도산대로               0.037         0.006                 1.3
4 20181231 동작대로               0.030         0.004                 0.9
5 20181231 시흥대로               0.056         0.006                 1.0
  아황산가스(ppm) 미세먼지(㎍/㎥) 초미세먼지(㎍/㎥)
1           0.005              52                31
2           0.005              50                29
3           0.004              48                27
4           0.008              48                26
5           0.005              44                28
                                                              주소
1 서울특별시 서초구 강남대로 201 서초구민회관 앞 중앙차로 (양재동)
2                         서울 성동구 강변북로 257 한강사업본부 옆
3                        서울 강남구 도산대로 104 신사역2번출구 앞
4                서울 동작구 동작대로 144 이수역 북단 버스중앙차로
5                     서울 금천구 독산동 996-9 한양수자인아파트 앞
                   지점명 sumdsum      lon      lat
1 강남대로(강남역-신분당)   80928 127.0360 37.48183
2                강변북로  230511 127.0416 37.53896
3        강남대로(신사역)   94804 127.0201 37.51605
4  동작대로(총신대입구역)  107203 126.9819 37.48649
5        시흥대로(시흥IC)   78455 126.8996 37.47557
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]

2.2.3 좌표계 정의

좌표계를 ’WGS1984’로 정의하기 위해 TA 지오코딩 세션에서 구득한 ’tl_emd_seoul_4326.shp’의 좌표계 정보를 사용하였습니다.

# OBTAIN CRS from existing shapefile
sigungu_shp <- readOGR(dsn = "./data", layer = "tl_emd._seoul_4326", encoding = "UTF-8", 
    stringsAsFactors = FALSE)

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Zulu-05\data", layer: "tl_emd._seoul_4326"
with 467 features
It has 3 fields

crs <- crs(sigungu_shp)

# define projection
coor <- SpatialPointsDataFrame(gc[, c(12:13)], gc, proj4string = crs)
tcoor <- SpatialPointsDataFrame(trafinfo1[, c(3:4)], trafinfo1, proj4string = crs)

# crs chk
crs(coor)

CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

crs(tcoor)

CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

잘 적용이 되었는지 plot을 해보았습니다.

# projection chk
plot(coor, main = "대기질관측소위치")

plot(tcoor, main = "교통량수집기위치")

2.2.4 Shape 저장 및 지도 도식

추후 사용 가능성을 염두에 두고 셰이프로 저장하고

# writeOGR(coor, paste0(getwd(),'/data/ex_shapefiles'), 'AQstations',
# driver='ESRI Shapefile', layer_options = 'encoding=UTF-8') writeOGR(tcoor,
# paste0(getwd(),'/data/ex_shapefiles'), 'traffic', driver='ESRI Shapefile',
# layer_options = 'encoding=UTF-8')

위에서 정제한 데이터로 ’tmap’라이브러리를 사용하여 도식하였습니다.

# echo map
tmap_mode("view")

tm_shape(coor) + tm_markers("측정소명", id = "측정소명", size = 0.5, shape = tmap_icons("./fig/lung.png"), 
    legend.col.show = FALSE) + tm_shape(tcoor) + tm_markers("지점명", id = "지점명", 
    size = 0.5, shape = tmap_icons("./fig/car.png"), legend.col.show = FALSE)

2.3.1 장소 선택

서울의 여러 관측소에서 나타난 일별 미세먼지 농도, 이산화질소 농도, 교통량 데이터 중 강변북로, 종로, 동작대로 관측소 데이터의 상관관계를 회귀분석을 통해 알아보기로 했습니다.

gangbyun <- subset(data_refined, 측정소명 == "강변북로")
jongno <- subset(data_refined, 측정소명 == "종로")
dongjak <- subset(data_refined, 측정소명 == "동작대로")

2.3.2 산점도 행렬 분석

강변북로, 종로, 동작대로의 여러 공기의 질 관련 데이터 중 산점도 행렬 분석(pairs.panels)을 이용해 각 변수와 교통량의 관계와 회귀분석 과정에서 나타나는 잔차의 분포에 대해서 알아봤습니다.

pairs.panels(gangbyun[, -c(1, 2, 9, 10)])

pairs.panels(jongno[, -c(1, 2, 9, 10)])

pairs.panels(dongjak[, -c(1, 2, 9, 10)])

잔차는 그래프를 보면 알 수 있듯이, 정확하지는 않지만 정규분포와 유사한 양상으로 나타났습니다.
따라서, 선형관계라는 가정이 적절하다는 것을 확인할 수 있습니다.

2.3.3 밀도함수

강변북로, 종로, 동작대로 세 장소 미세먼지 농도의 측정값을 바탕으로 밀도함수를 그려봤습니다.

ggplot(jongno, aes(`미세먼지(㎍/㎥)`)) + geom_density(kernel = "gaussian")

ggplot(gangbyun, aes(`미세먼지(㎍/㎥)`)) + geom_density(kernel = "gaussian")

ggplot(dongjak, aes(`미세먼지(㎍/㎥)`)) + geom_density(kernel = "gaussian")

2.3.4 미세먼지-교통량 회귀분석

회귀분석을 통해 강변북로, 종로, 강남대로 세 장소의 교통량과 미세먼지 농도의 상관관계에 대해 알아봤습니다.

gf_point(gangbyun$`미세먼지(㎍/㎥)` ~ gangbyun$sumdsum, data = gangbyun) %>% gf_lm(interval = "prediction", 
    fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 215000, 
    y = 110, label = "y = 0.0008711x -167.4") + geom_text(x = 215000, y = 102, label = "p-value = 7.707e-09, R²= 0.1064") + 
    labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")

gf_point(jongno$`미세먼지(㎍/㎥)` ~ jongno$sumdsum, data = jongno) %>% gf_lm(interval = "prediction", 
    fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 35000, 
    y = 100, label = "y = 0.001892x -51.91") + geom_text(x = 35000, y = 92, label = "p-value = 0.0001899, R²= 0.09567") + 
    labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")

gf_point(dongjak$`미세먼지(㎍/㎥)` ~ dongjak$sumdsum, data = dongjak) %>% gf_lm(interval = "prediction", 
    fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 1e+05, 
    y = 105, label = "y = 0.001169x -85.21") + geom_text(x = 1e+05, y = 97, label = "p-value = 0.0002134, R²= 0.08916") + 
    labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")

2.3.5 미세먼지-교통량 회귀분석 검증

lm(gangbyun$`미세먼지(㎍/㎥)` ~ gangbyun$sumdsum) %>% summary()

Call:
lm(formula = gangbyun$`미세먼지(㎍/㎥)` ~ gangbyun$sumdsum)

Residuals:
   Min     1Q Median     3Q    Max 
-46.84 -17.64  -3.32  15.22  77.22 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -1.674e+02  3.557e+01  -4.707 3.85e-06 ***
gangbyun$sumdsum  8.711e-04  1.465e-04   5.946 7.71e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23.21 on 297 degrees of freedom
Multiple R-squared:  0.1064,    Adjusted R-squared:  0.1034 
F-statistic: 35.35 on 1 and 297 DF,  p-value: 7.707e-09

lm(jongno$`미세먼지(㎍/㎥)` ~ jongno$sumdsum) %>% summary()

Call:
lm(formula = jongno$`미세먼지(㎍/㎥)` ~ jongno$sumdsum)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.886 -11.990  -2.771   8.452  81.063 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -5.191e+01  2.227e+01  -2.330  0.02123 *  
jongno$sumdsum  1.892e-03  4.934e-04   3.835  0.00019 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.24 on 139 degrees of freedom
Multiple R-squared:  0.09567,   Adjusted R-squared:  0.08916 
F-statistic:  14.7 on 1 and 139 DF,  p-value: 0.0001899

lm(dongjak$`미세먼지(㎍/㎥)` ~ dongjak$sumdsum) %>% summary()

Call:
lm(formula = dongjak$`미세먼지(㎍/㎥)` ~ dongjak$sumdsum)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.061 -17.583  -3.873  11.453  90.802 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -8.521e+01  3.501e+01  -2.434 0.015557 *  
dongjak$sumdsum  1.169e-03  3.117e-04   3.751 0.000213 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 23.86 on 284 degrees of freedom
Multiple R-squared:  0.04721,   Adjusted R-squared:  0.04385 
F-statistic: 14.07 on 1 and 284 DF,  p-value: 0.0002134

gf_point(gangbyun$`이산화질소농도(ppm)` ~ gangbyun$sumdsum, data = gangbyun) %>% 
    gf_lm(interval = "prediction", fill = "red") %>% gf_lm(interval = "confidence", 
    fill = "navy") + geom_text(x = 215000, y = 0.07, label = "y = 0.00000055789x -0.09588") + 
    geom_text(x = 215000, y = 0.065, label = "p-value = 7.178e-10, R²= 0.1172") + 
    labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")

gf_point(jongno$`이산화질소농도(ppm)` ~ jongno$sumdsum, data = jongno) %>% gf_lm(interval = "prediction", 
    fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 35000, 
    y = 0.065, label = "y = 0.000000869x -0.002781") + geom_text(x = 35000, y = 0.06, 
    label = "p-value = 0.01181, R²=0.03786") + labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")

gf_point(dongjak$`이산화질소농도(ppm)` ~ dongjak$sumdsum, data = dongjak) %>% gf_lm(interval = "prediction", 
    fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 1e+05, 
    y = 0.08, label = "y = 0.0000008714x -0.06077") + geom_text(x = 1e+05, y = 0.076, 
    label = "p-value = 1.463e-05, R²= 0.06082") + labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")

2.3.6 이산화질소-교통량 회귀분석 검증

lm(gangbyun$`이산화질소농도(ppm)` ~ gangbyun$sumdsum) %>% summary()

Call:
lm(formula = gangbyun$`이산화질소농도(ppm)` ~ gangbyun$sumdsum)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.037627 -0.008847  0.001063  0.009473  0.038892 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -9.588e-02  2.126e-02   -4.51 9.33e-06 ***
gangbyun$sumdsum  5.578e-07  8.756e-08    6.37 7.18e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01387 on 297 degrees of freedom
Multiple R-squared:  0.1202,    Adjusted R-squared:  0.1172 
F-statistic: 40.58 on 1 and 297 DF,  p-value: 7.178e-10

lm(jongno$`이산화질소농도(ppm)` ~ jongno$sumdsum) %>% summary()

Call:
lm(formula = jongno$`이산화질소농도(ppm)` ~ jongno$sumdsum)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.028371 -0.009724 -0.000216  0.008137  0.032833 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)    -2.781e-03  1.538e-02  -0.181   0.8568  
jongno$sumdsum  8.690e-07  3.406e-07   2.551   0.0118 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01259 on 139 degrees of freedom
Multiple R-squared:  0.04473,   Adjusted R-squared:  0.03786 
F-statistic: 6.509 on 1 and 139 DF,  p-value: 0.01181

lm(dongjak$`이산화질소농도(ppm)` ~ dongjak$sumdsum) %>% summary()

Call:
lm(formula = dongjak$`이산화질소농도(ppm)` ~ dongjak$sumdsum)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.02460 -0.01224 -0.00399  0.01151  0.03884 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -6.077e-02  2.219e-02  -2.739  0.00656 ** 
dongjak$sumdsum  8.714e-07  1.976e-07   4.411 1.46e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01512 on 284 degrees of freedom
Multiple R-squared:  0.06412,   Adjusted R-squared:  0.06082 
F-statistic: 19.46 on 1 and 284 DF,  p-value: 1.463e-05

3.2.1 미세먼지

우선, 2018년 데이터 중 일요일의 대기질 데이터를 추출하였습니다.

# Extract Sunday air quality data
start_date <- as.Date("2018-01-1")
end_date <- as.Date("2018-12-31")
df <- data.frame(seq(as.Date(start_date), as.Date(end_date), by = 1))
names(df) <- "date_18"
df$weekday <- weekdays(as.Date(df$date_18))

str(df)

'data.frame':   365 obs. of  2 variables:
 $ date_18: Date, format: "2018-01-01" "2018-01-02" ...
 $ weekday: chr  "Monday" "Tuesday" "Wednesday" "Thursday" ...

sunday <- filter(df, weekday == "Sunday")

datelist <- air$측정일시
date <- ymd(datelist)
air <- air %>% mutate(date)

sunday_air <- air %>% filter(date %in% sunday$date_18)
sunday_air <- rename(sunday_air, 측정일시 = "측정일시", 미세먼지 = "미세먼지(㎍/㎥)", 
    초미세먼지 = "초미세먼지(㎍/㎥)", 이산화질소 = "이산화질소농도(ppm)", 오존 = "오존농도(ppm)", 
    이산화탄소 = "이산화탄소농도(ppm)", 아황산가스 = "아황산가스(ppm)")

head(sunday_air)

# A tibble: 6 x 9
  측정일시 측정소명 이산화질소  오존 이산화탄소 아황산가스 미세먼지 초미세먼지
     <dbl> <chr>         <dbl> <dbl>      <dbl>      <dbl>    <dbl>      <dbl>
1 20180107 강남구        0.048 0.003        0.7      0.006       46         42
2 20180107 강남대로      0.065 0.003        1.1      0.006       69         31
3 20180107 강동구        0.049 0.004        0.9      0.005       57         45
4 20180107 강변북로      0.041 0.004        0.9      0.005       54         25
5 20180107 강북구        0.047 0.007        0.9      0.004       49         31
6 20180107 강서구        0.056 0.004        1        0.004       67         45
# ... with 1 more variable: date <date>

일요일 대기질 데이터를 측정일시를 기준으로 그룹화하여 결측치를 제거한 후, 미세먼지 농도의 일별 평균을 구했습니다.

# mean of fine dust
n <- sunday_air %>% group_by(측정일시) %>% summarise(n = n())

mean <- sunday_air %>% group_by(측정일시) %>% summarise(미세먼지평균 = mean(미세먼지, 
    na.rm = TRUE))
mean_table <- inner_join(mean, n, by = "측정일시")

summary(mean_table$미세먼지평균)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  12.33   23.02   34.00   37.73   52.00   77.36 

초미세먼지 농도의 경우도 마찬가지로 평균을 구한 뒤 합쳤습니다.

## mean of ultra fine dust
mean2 <- sunday_air %>% group_by(측정일시) %>% summarise(초미세먼지평균 = mean(초미세먼지, 
    na.rm = TRUE))
mean_table <- inner_join(mean2, mean_table, by = "측정일시")

summary(mean_table$초미세먼지평균)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  5.872  11.128  20.308  21.647  29.000  50.538 

head(mean_table)

# A tibble: 6 x 4
  측정일시 초미세먼지평균 미세먼지평균     n
     <dbl>          <dbl>        <dbl> <int>
1 20180107           37.9         57.5    39
2 20180114           50.5         71.1    39
3 20180121           34.0         65.9    39
4 20180128           19.2         34      39
5 20180204           21           33.6    39
6 20180211           14.7         53.3    39

summary(mean_table)

    측정일시        초미세먼지평균    미세먼지평균         n        
 Min.   :20180107   Min.   : 5.872   Min.   :12.33   Min.   :39.00  
 1st Qu.:20180415   1st Qu.:11.128   1st Qu.:23.02   1st Qu.:39.00  
 Median :20180729   Median :20.308   Median :34.00   Median :39.00  
 Mean   :20180703   Mean   :21.647   Mean   :37.73   Mean   :40.31  
 3rd Qu.:20181014   3rd Qu.:29.000   3rd Qu.:52.00   3rd Qu.:39.00  
 Max.   :20181230   Max.   :50.538   Max.   :77.36   Max.   :46.00  

3.2.2 시청률

시청률 데이터 역시 결측치를 제거한 후 조인합니다.

# viewrate data
viewrate <- read_csv("./data/viewratedata.csv")
viewrate <- viewrate %>% mutate(total = as.numeric(viewrate$Total)) %>% dplyr::select(Date, 
    total)

# eliminate missing values
totalrate <- rename(viewrate, 측정일시 = "Date", total_viewrate = "total")
totalrate <- na.omit(totalrate)


head(totalrate)

# A tibble: 6 x 2
  측정일시 total_viewrate
     <dbl>          <dbl>
1 20180107           36.3
2 20180114           35.5
3 20180121           31.4
4 20180128           35.2
5 20180204           32.8
6 20180218           33.1

summary(totalrate)

    측정일시        total_viewrate
 Min.   :20180107   Min.   :21.1  
 1st Qu.:20180418   1st Qu.:26.0  
 Median :20180708   Median :28.0  
 Mean   :20199832   Mean   :28.3  
 3rd Qu.:20181010   3rd Qu.:30.8  
 Max.   :21080325   Max.   :36.3  

3.2.3 완성된 데이터

준비 완료!

# data join
mean_table <- inner_join(mean_table, totalrate, by = "측정일시")
summary(mean_table)

    측정일시        초미세먼지평균    미세먼지평균         n        
 Min.   :20180107   Min.   : 5.872   Min.   :12.33   Min.   :39.00  
 1st Qu.:20180420   1st Qu.:11.215   1st Qu.:23.24   1st Qu.:39.00  
 Median :20180726   Median :21.346   Median :35.01   Median :39.00  
 Mean   :20180713   Mean   :22.395   Mean   :38.46   Mean   :40.48  
 3rd Qu.:20181016   3rd Qu.:29.530   3rd Qu.:52.17   3rd Qu.:39.00  
 Max.   :20181230   Max.   :50.538   Max.   :77.36   Max.   :46.00  
 total_viewrate 
 Min.   :21.10  
 1st Qu.:26.27  
 Median :28.35  
 Mean   :28.61  
 3rd Qu.:31.23  
 Max.   :36.30  

DescTools::Desc(mean_table)

------------------------------------------------------------------------------ 
Describe mean_table (tbl_df, tbl, data.frame):

data frame: 40 obs. of  5 variables
        40 complete cases (100.0%)

  Nr  ColName         Class    NAs  Levels
  1   측정일시        numeric  .          
  2   초미세먼지평균  numeric  .          
  3   미세먼지평균    numeric  .          
  4   n               integer  .          
  5   total_viewrate  numeric  .          


------------------------------------------------------------------------------ 
1 - 측정일시 (numeric)

    length         n       NAs    unique        0s      mean    meanCI
        40        40         0       = n         0  2.02e+07  2.02e+07
              100.0%      0.0%                0.0%            2.02e+07
                                                                      
       .05       .10       .25    median       .75       .90       .95
  2.02e+07  2.02e+07  2.02e+07  2.02e+07  2.02e+07  2.02e+07  2.02e+07
                                                                      
     range        sd     vcoef       mad       IQR      skew      kurt
  1'123.00    366.00      0.00    444.78    595.50     -0.19     -1.34
                                                                      
lowest : 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07
highest: 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07

------------------------------------------------------------------------------ 
2 - 초미세먼지평균 (numeric)

     length          n        NAs     unique         0s       mean     meanCI
         40         40          0        = n          0  22.395112  18.566333
                100.0%       0.0%                  0.0%             26.223891
                                                                             
        .05        .10        .25     median        .75        .90        .95
   6.787821   9.615824  11.215476  21.346154  29.529762  39.550343  44.500328
                                                                             
      range         sd      vcoef        mad        IQR       skew       kurt
  44.666667  11.971837   0.534574  14.660090  18.314286   0.591248  -0.576656
                                                                             
lowest : 5.871795, 6.166667, 6.820513, 7.542857, 9.846154
highest: 39.456522, 40.394737, 44.325581, 47.820513, 50.538462

------------------------------------------------------------------------------ 
3 - 미세먼지평균 (numeric)

    length         n       NAs    unique        0s      mean    meanCI
        40        40         0       = n         0  38.45879  32.55534
              100.0%      0.0%                0.0%            44.36224
                                                                      
       .05       .10       .25    median       .75       .90       .95
  13.06855  15.46923  23.24240  35.01429  52.17391  66.16410  69.10513
                                                                      
     range        sd     vcoef       mad       IQR      skew      kurt
  65.02564  18.45892   0.47997  18.91139  28.93152   0.41535  -0.98042
                                                                      
lowest : 12.33333, 12.69444, 13.08824, 14.38462, 15.58974
highest: 65.92308, 68.33333, 69.0, 71.10256, 77.35897

------------------------------------------------------------------------------ 
4 - n (integer)

  length       n    NAs  unique     0s   mean  meanCI
      40      40      0       3      0  40.48   39.59
          100.0%   0.0%           0.0%          41.36
                                                     
     .05     .10    .25  median    .75    .90     .95
   39.00   39.00  39.00   39.00  39.00  46.00   46.00
                                                     
   range      sd  vcoef     mad    IQR   skew    kurt
    7.00    2.78   0.07    0.00   0.00   1.29   -0.30
                                                     

   level  freq   perc  cumfreq  cumperc
1     39    31  77.5%       31    77.5%
2     45     4  10.0%       35    87.5%
3     46     5  12.5%       40   100.0%

heap(?): remarkable frequency (77.5%) for the mode(s) (= 39)

------------------------------------------------------------------------------ 
5 - total_viewrate (numeric)

  length       n     NAs  unique      0s    mean  meanCI
      40      40       0      39       0  28.613  27.458
          100.0%    0.0%            0.0%          29.767
                                                        
     .05     .10     .25  median     .75     .90     .95
  23.680  24.810  26.275  28.350  31.225  33.200  35.215
                                                        
   range      sd   vcoef     mad     IQR    skew    kurt
  15.200   3.610   0.126   3.781   4.950   0.153  -0.447
                                                        
lowest : 21.1, 21.4, 23.8, 24.0, 24.9
highest: 33.1, 34.1, 35.2, 35.5, 36.3

DT::datatable(mean_table, rownames = F, caption = "")

3.3.1 변수 검정

분석에 앞서 변수들의 정규성을 검정하였고, 검정 결과 미세먼지와 초미세먼지의 p-value가 유의 수준을 넘지 않았으므로 귀무가설 기각에 실패하였습니다.

# shapiro test
shapiro.test(mean_table$미세먼지평균)  #p-value 0.05691

    Shapiro-Wilk normality test

data:  mean_table$미세먼지평균
W = 0.94637, p-value = 0.05691

shapiro.test(mean_table$초미세먼지평균)  #p-value 0.04389

    Shapiro-Wilk normality test

data:  mean_table$초미세먼지평균
W = 0.94306, p-value = 0.04389

shapiro.test(mean_table$total_viewrate)  #p-value 0.8075

    Shapiro-Wilk normality test

data:  mean_table$total_viewrate
W = 0.98326, p-value = 0.8075

따라서 대기질 데이터에 log를 취하여 정규성을 갖도록 하였습니다.

# normalize variables
mean_table <- transform(mean_table, fine_log = log(미세먼지평균), ultrafine_log = log(초미세먼지평균))

그 결과는 다음과 같습니다.

# again shapiro test
shapiro.test(mean_table$fine_log)  #p-value= 0.1346

    Shapiro-Wilk normality test

data:  mean_table$fine_log
W = 0.95723, p-value = 0.1346

shapiro.test(mean_table$ultrafine_log)  #p-value = 0.2193

    Shapiro-Wilk normality test

data:  mean_table$ultrafine_log
W = 0.96345, p-value = 0.2193

3.3.2 가설설정

설정한 가설은
H0 : 미세먼지 농도는 시청률에 영향을 주지 않는다.
H1 : 미세먼지 농도는 시청률에 영향을 준다.

3.3.3 상관분석

회귀식을 도출하기에 앞서, cortest를 진행하였습니다. 그 결과 상관계수는 0.3541이고, p-value는 0.02494이기 때문에 두 변수간 유의한 상관관계가 있음을 알 수 있습니다.

# correlation

cor.test(mean_table$fine_log, mean_table$total_viewrate)

    Pearson's product-moment correlation

data:  x and y
t = 2.3348, df = 38, p-value = 0.02494
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.04797772 0.59955155
sample estimates:
      cor 
0.3541935 

3.3.4 단일 회귀 분석

미세먼지와 시청률간의 회귀분석을 진행하였습니다.

# regression of fine dust
model1 <- lm(total_viewrate ~ fine_log, data = mean_table)

summary(model1)

Call:
lm(formula = total_viewrate ~ fine_log, data = mean_table)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.8905 -1.7322 -0.4328  1.6073  6.5856 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   20.000      3.728   5.364 4.24e-06 ***
fine_log       2.443      1.046   2.335   0.0249 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.42 on 38 degrees of freedom
Multiple R-squared:  0.1255,    Adjusted R-squared:  0.1024 
F-statistic: 5.451 on 1 and 38 DF,  p-value: 0.02494

회귀 분석을 진행한 결과,
1. F-Statistic의 p-value 값이 유의 수준인 0.05보다 작게 도출되었습니다.

2. Coefficients p-value값이 0.024로 이 역시 유의 수준 0.05보다 작게 도출되었습니다.
   

p-value에 따르면 미세먼지 농도는 시청률에 유의한 영향을 주는 것으로 해석할 수 있습니다.

3. 그러나 결정계수의 값이 0.1수준으로 낮게 나왔기 때문에 이 회귀 모형은 통계적으로 완벽히 유의하다고 할 수 없습니다.

stargazer(model1, title = "total_ viewrate ~ fine_dust", align = TRUE, type = "html")

gf_point(total_viewrate ~ fine_log, data = mean_table) %>% gf_lm(interval = "confidence") + 
    xlab("fine dust(㎍/㎥)") + ylab("total viewrate(%)")

혹시나 하는 마음에 초미세먼지 데이터로도 회귀분석을 진행하였습니다.

# regression of ultrafinedust
model2 <- lm(total_viewrate ~ ultrafine_log, data = mean_table)

summary(model2)

Call:
lm(formula = total_viewrate ~ ultrafine_log, data = mean_table)

Residuals:
   Min     1Q Median     3Q    Max 
-5.766 -1.939 -0.319  1.860  6.590 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    20.3200     2.6966   7.535 4.69e-09 ***
ultrafine_log   2.8053     0.8954   3.133  0.00333 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.26 on 38 degrees of freedom
Multiple R-squared:  0.2053,    Adjusted R-squared:  0.1844 
F-statistic: 9.815 on 1 and 38 DF,  p-value: 0.003327

역시나 비슷한 결과가 도출되었습니다.

gf_point(total_viewrate ~ ultrafine_log, data = mean_table) %>% gf_lm(interval = "confidence") + 
    xlab("ultrafine dust(㎍/㎥)") + ylab("total viewrate(%)")

3.3.5 원인 찾기

plot을 통해 회귀 분석 결과를 자세히 보겠습니다.

plot(model1)

plot 상으로는 큰 문제가 없어 보입니다. 따라서 보다 정밀하게 알아보고자

gvlma패키지를 통해 다시 살펴보면, 역시 Heteroscedasticity 검정에서 통과하지 못했음을 확인할 수 있습니다.

summary(gvlma(model1))  #omg....

Call:
lm(formula = total_viewrate ~ fine_log, data = mean_table)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.8905 -1.7322 -0.4328  1.6073  6.5856 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   20.000      3.728   5.364 4.24e-06 ***
fine_log       2.443      1.046   2.335   0.0249 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.42 on 38 degrees of freedom
Multiple R-squared:  0.1255,    Adjusted R-squared:  0.1024 
F-statistic: 5.451 on 1 and 38 DF,  p-value: 0.02494


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = model1) 

                      Value p-value                   Decision
Global Stat        7.763782 0.10062    Assumptions acceptable.
Skewness           0.002025 0.96410    Assumptions acceptable.
Kurtosis           0.262686 0.60828    Assumptions acceptable.
Link Function      0.115021 0.73450    Assumptions acceptable.
Heteroscedasticity 7.384050 0.00658 Assumptions NOT satisfied!

plot을 통해 확인한 outlier를 제거한 후 다시 회귀분석을 해보았습니다.

# remove outlier
fit_table <- mean_table[-c(1, 4, 10, 16, 17, 27), ]
re <- lm(total_viewrate ~ log(미세먼지평균), data = fit_table)

# again regression
plot(re)

g1 <- gvlma(re)
summary(g1)

Call:
lm(formula = total_viewrate ~ log(미세먼지평균), data = fit_table)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2213 -1.5737 -0.2117  1.1307  4.1737 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         15.783      3.055   5.166 1.23e-05 ***
log(미세먼지평균)    3.653      0.865   4.223 0.000186 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.49 on 32 degrees of freedom
Multiple R-squared:  0.3579,    Adjusted R-squared:  0.3378 
F-statistic: 17.84 on 1 and 32 DF,  p-value: 0.0001863


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = re) 

                     Value p-value                Decision
Global Stat        3.78793 0.43547 Assumptions acceptable.
Skewness           0.04969 0.82360 Assumptions acceptable.
Kurtosis           0.33388 0.56338 Assumptions acceptable.
Link Function      0.01352 0.90745 Assumptions acceptable.
Heteroscedasticity 3.39084 0.06556 Assumptions acceptable.

summary(re)

Call:
lm(formula = total_viewrate ~ log(미세먼지평균), data = fit_table)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.2213 -1.5737 -0.2117  1.1307  4.1737 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         15.783      3.055   5.166 1.23e-05 ***
log(미세먼지평균)    3.653      0.865   4.223 0.000186 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.49 on 32 degrees of freedom
Multiple R-squared:  0.3579,    Adjusted R-squared:  0.3378 
F-statistic: 17.84 on 1 and 32 DF,  p-value: 0.0001863

그 결과 p-value값은 더욱 작아졌고, 위에서 문제가 되었던 R Squared 값이 0.3378정도로 증가하였음을 확인 할 수 있었습니다. 설명변수로 반응변수를 약 33% 설명할 수 있습니다.

4.2.1 가설 설정

설정한 가설은 다음과 같습니다
HO : 강수량은 미세먼지 농도에 영향을 주지 않는다
H1 : 강수량은 미세먼지 농도에 영향을 준다

4.2.2 회귀분석

왼쪽으로 치우쳐진 분포를 보이는 일강수량과 미세먼지를 log 변환하였습니다.

# normalize variables
rain$lograin = with(rain, log(일강수량))
rain$logdust = with(rain, log(미세먼지))

par(mfrow = c(1, 2))
densityplot(~logdust, data = rain)

densityplot(~lograin, data = rain)

이후 회귀 분석을 시행하였고, 회귀 분석의 결과는 다음과 같습니다.

gf_point(logdust ~ lograin, data = rain)

gf_point(logdust ~ lograin, data = rain) %>% gf_lm(interval = "prediction") %>% gf_lm(interval = "confidence")

rain$lograin[is.infinite(rain$lograin)] <- NA
lm_rain = lm(logdust ~ lograin, data = rain)
summary(lm_rain)

Call:
lm(formula = logdust ~ lograin, data = rain)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.23179 -0.44525  0.01916  0.39422  1.31442 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.32133    0.07663  43.340  < 2e-16 ***
lograin     -0.13392    0.03451  -3.881  0.00021 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5811 on 82 degrees of freedom
  (21 observations deleted due to missingness)
Multiple R-squared:  0.1552,    Adjusted R-squared:  0.1449 
F-statistic: 15.06 on 1 and 82 DF,  p-value: 0.0002096

p-value 값이 0.05보다 작은 0.0002 정도이므로 회귀분석 결과 귀무가설을 기각 합니다. 즉, 일 강수량은 미세먼지 농도에 영향을 준다고 볼 수 있습니다. 그러나, 결정계수 R-squared값이 약 0.14 정도로 다소 낮다는 한계를 가집니다.